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Abstract 

The transport and gain properties of quantum cascade (QC) structures are investigated using 
a nonequilibrium Green's function (NGF) theory which includes quantum effects beyond a Boltz- 
mann transport description. In the NGF theory, we include interface roughness, impurity, and 
electron-phonon scattering processes within a self-consistent Born approximation, and electron- 
electron scattering in a mean-field approximation. With this theory we obtain a description of the 
nonequilibrium stationary state of QC structures under an applied bias, and hence we determine 
transport properties, such as the current-voltage characteristic of these structures. We define two 
contributions to the current, one contribution driven by the scattering- free part of the Hamilto- 
nian, and the other driven by the scattering Hamiltonian. We find that the dominant part of the 
current in these structures, in contrast to simple superlattice structures, is governed mainly by 
the scattering Hamiltonian. In addition, by considering the linear response of the stationary state 
of the structure to an applied optical field, we determine the linear susceptibility, and hence the 
gain or absorption spectra of the structure. A comparison of the spectra obtained from the more 
rigorous NGF theory with simpler models shows that the spectra tend to be offset to higher values 
in the simpler theories. 
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I. INTRODUCTION 



Quantum cascade (QC) structures are semiconductor heterostructures grown with a com- 
plicated sequence of alternating layers of different semiconductor materials and with varying 
thicknesses. This sequence of layers is repeated many times, up to tens or even over a hun- 
dred periods. Figure [I] shows an example of the conduction-band line-up in a QC structure. 
These structures form the basis of a new type of semiconductor laser ,@ in which the laser light 
emission occurs through intersubband or interminiband transitions (in most cases within the 
conduction band) rather than interband transitions. These lasers have a great variety of de- 
signs, and a recent review is given in Ref. 0. Until recently^ all quantum cascade laser (QCL) 
structures were designed so that each period in the structure contains an active region in 
which the lasing transition occurs, and a separate injector region. The injector acts as a 
reservoir of electrons for injection into the active region of the next stage. It also acts as a 
collector of electrons from the preceding active region. The direction of the electron flow is 
indicated in Fig. [I], and the electron flow is seen to resemble a cascade as the electrons move 
from one stage to the next when a bias is applied. Hence, the electron transport through 
a QCL structure is a complicated interplay between relaxation through light emission and 
scattering in the active region, and transmission through tunneling and scattering in the 
injector region.i 

Initially, most theoretical investigations&i! of transport in QC structures have focused on 
the role played by scattering processes in determining transport properties and the dynamics 
of the electron distributions in these structures. Very recently, however, the first theoretical 
investigations of quantum transport in these structures which have considered or incor- 
porated quantum effects beyond a Boltzmann equation approach have been reported, 
theoretical study of quantum transport may be treated using the density matrix formalismJ§ 
or with a nonequilibrium Green's function (NGF) approach£2Elll'El 

In the work reported here, we extend the NGF theory described in Ref. [T^ to the study 
of quantum transport in QC structures. Very early results from this investigation have been 
reported in Refs. |7| and [T3|. In this paper, we present further and more detailed results from 
this study. In addition, we have extended the work further to the problem of evaluating the 
gain or absorption spectra of these structures,^ and we also report and discuss results of 
this work here. 

In the following section, we describe the theoretical formulation that we use to derive the 
transport properties, and then the linear optical response of QC structures. In Sec. |TTT|, we 
apply this theory to example QC structures, and describe the results obtained. The last 
section contains a summary and conclusion. 



II. THEORETICAL FORMULATION 
A. Basis states and Hamiltonian 

We model the QC structure as a periodic superlattice structure, in which each period 
contains N s semiconductor layers with varying thicknesses. The Hamiltonian H which we 
use to model this system may be separated into two parts: H = H Q + H SC a,tt- H Q contains 
the superlattice potential and a static electric field £ applied in the growth direction, i.e., 
H = Hsl + H £ . The Hamiltonian H scat t describes the scattering processes included in the 
theory. 
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The Hamiltonian is expressed in a set of basis states, \]/k,a( r ; z) — (e lk ' r / VA)ip a (z) , which 
we assume separable, although this is an approximation when the effective mass is position 
dependent. In the plane of the semiconductor layers, the basis functions behave as plane 
waves. The normalization constant A is the sample area in this plane. The position vector r 
and in-plane wave vector k are two-dimensional (2D) vectors. In the growth direction, here 
labeled z, there are several possible choices (see Ref. [13] for a discussion) for the functions 
ip a {z)- Although the physical results obtained from the theory should be independent of the 
specific choice of these functions, there are different advantages or drawbacks attached to a 
given choice. For instance, as we discuss next, one choice may be more suited to expedit- 
ing the numerical computation, while another choice may more easily allow the extraction 
of physical information in a form that can be compared with experimental measurements. 
Possible choices are (i) Bloch functions which are eigenstates of the bare superlattice poten- 
tial Hsl, and are spatially extended across the whole structure, (ii) Wannier-Stark states 
which are eigenstates of fL, i.e., of the superlattice potential and the applied bias. These 
eigenstates are metastableE3 (non-Hermitian) with complex energies. They are often treated 
approximately as stationary states (their metastable nature is neglected), and this leads to 
an ambiguity in the definition of these states depending on how this approximation is made, 
(iii) Wannier functions, which should not be confused with Wannier-Stark functions, may be 
constructed so that they are spatially well localized. In particular, they may be constructed 
as eigenstates of the position operator. In symmetric superlattice structures, these eigen- 
states are maximally localized.li3 The spatial localization of the Wannier states enables one 
to setup a picture of transport in which scattering occurs from one spatial region of the struc- 
ture to another. The Wannier functions do not depend on bias, and this is computationally 
advantageous since coupling matrix elements between Wannier functions may be calculated 
only once at zero bias, and the same matrix elements may then be used at any applied bias. 
For these reasons, the theory presented here is formulated in the Wannier function basis 
(see Appendix [A]), which we derive as eigenstates of the position operator. A disadvantage 
in using this basis is that the Wannier functions are not energy eigenstates of the struc- 
ture, and it is therefore difficult to obtain a physical interpretation in the energy picture or 
domain. For instance, there is no one-to-one correspondence between an optical transition 
in the structure and a transition between a pair of Wannier functions. This difficulty may 
be partially circumvented, however, at a later stage in the calculation, by transforming the 
results obtained into the Wannier-Stark basis. We emphasize again that the basis choice is 
not in itself an important issue, i.e., the physical content of the theory should not depend 
on this choice, but a suitable choice of basis can facilitate the numerical computation (e.g., 
Wannier states), or more easily allow physical interpretation (e.g., Wannier-Stark states). 

Expressing the Hamiltonian in the Wannier basis we obtain 

n,u k,s 

(^n+l,k,s^n,k,s + ^n-l,k,s^n,k, J] > (1) 

where the index n labels a period in the superlattice, and the index v labels a Wannier 
function ip u (z) within a period. c^f ks and a^ ks are creation and annihilation operators for 
an electron with in-plane wave vector k, and spin index s, in the vth Wannier level, in 
period n. As stated earlier, the Wannier functions are not eigenstates of Hsl, and hence T" 
represents the off-diagonal couplings between Wannier levels in different periods, and 
represents the diagonal elements of Hsl in this basis. We keep only terms in T%, i.e., we 
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consider only couplings between adjacent periods. The next-nearest-neighbor couplings Tg 
are two or more orders of magnitude smaller. The Hamiltonian H^, due to the electric field 
S, is written as 

— eSRi l&n+i^s «n,k,s + ^n|k,s ^n+l,k,s ]}) (2) 

where i?f = f dz ip^(z — Id) zip u (z). d is the length of one period and e < is the electron 
charge. 

In the scattering Hamiltonian H sca ,tt, we include interface roughness, impurity, and 
electron-phonon scattering processes. Both acoustic phonons and longitudinal optical (LO) 
phonons are considered. The electron-acoustic-phonon scattering rates are far smaller than 
those of the other scattering processes (Table |), but we include these to provide a channel 
for small energy transfers between the electrons and lattice. This opens up a route for the 
carrier-lattice system to move towards thermal equilibrium, particularly at low bias when 
the electrons may have insufficient energy to emit LO-phonons. 

In the NGF theory, the scattering processes expressed in H scatt are treated in the form of 
self-energies which are described below. We have also included electron-electron scattering 
in a mean-field approximation so that the interaction between electrons appears as an ad- 
ditional single-particle potential in H scatt , i.e., we replace F scatt with # s ' ca tt = #scatt + H MF . 
Hmf is written as 

^MF = 51 X^ MF W>™' «m,k,s ^n,k,s> (3) 
mn,nv k.s 

where [T4iF]m/j,w = / dz ip^Jz) Vmf( z ) ipnui 2 )- The evaluation of the mean-field potential 
Vmf(z) is described in the following subsection. 



B. Quantum transport equations and self-energies 

To determine the transport properties of a QC structure under applied bias, we need first 
to obtain a description of the nonequilibrium stationary state of the system. This information 
is contained in the nonequilibrium Green's functions: the retarded Green's function G iet (E) 
and the correlation function G < (E) (bold type G, and £ below, indicates a matrix in, e.g., 
the Wannier basis). The quantum transport equations, which these functions obey, can be 
derived00 from the Hamiltonian H and have the following form. The Dyson equation for 
the matrix element G^ tt k(-^) * s written as 

EG aia 2 ,k( E ) - [( H o + #MF)ai/3,k + ^0,k( E ) G % 2 ,k( E ) = <W 2 > ( 4 ) 
P 

where ati, «2, and (3 are general indices that include both the period and Wannier level 
indices, e.g. cm = (m,fi). The correlation function G^ lQ , 2 ^(E) is obtained from the Keldysh 
relation: 

G a ia2 ,k( E ) = G ai0,k( E ) ^Pf3>,k( E ) G f3'Z 2 ,k( E )i ( 5 ) 

where G adv (E) = [G ret (E)}l 
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The self-energies S rct (i?) and 'S < (E) originate from the scattering processes contained in 
H sca ,tt (i-e., excluding Huf), and are evaluated within the self-consistent Born approximation. 
Assuming that the diagonal parts of the Green's functions and the self-energies dominate 
(because the basis states are spatially localized), the self-energy for interface roughness or 
impurity scattering has the form 

^::f imp (E) = E(lC; sVimP (k - k')D G<^{E). (6) 

V™p S (k' — k) and ^"^(k' — k) are matrix elements for interface roughness and impurity 
scattering. In the former case, the angle brackets denote an average over all possible distribu- 
tions of variations in the interface thickness, and in the latter case these brackets denote an 
average over all possible impurity distributions. The equation for $] ret ' rough / imp (£') j s obtained 
by making the replacements: s < ' rou s h / imp ( J B) -> s ret ' rough / imp (£;), and G < (E) -> G TCt (E). 
For optical or acoustic phonon scattering, the self-energies are 

= E l<7( k > k ')l 2 

j9,k' 



X 



[/fl(-E'phoii) + 1] ^/3/3,k' (-^ — -Ephon 



+ /s(-Ephon) G^ ik ,(E + -Ephon) 



2 

+ % 



and 



/3,k' 



x j/B^phon) Gpf3,k'{E — Sphou) 

+ [/e(-Ephon) + 1] Gppy^E + S p hon)|, (8) 



where V phon represents the interaction with optical or acoustic phonons, and E phon rep- 
resents the energy of the optical or acoustic phonon. fs(E) = l/fexp^/fcgT) — 1] is the 
equilibrium phonon distribution at a temperature T. This is the only place where the lattice 
temperature appears in this theory. Further detail about the evaluation of the scattering 
matrix elements and the self-energies is given in Appendix |B[ For simplicity in the nu- 
merical evaluation, we neglect the last line containing the principle value terms in Eq. (^). 
To further expedite the numerical computation, we use momentum-independent scattering 
matrix elements V^ on / rough / imp (k typ , k' typ ), employing a representative momentum k typ (see 
Appendix ||). 

As mentioned in the previous section, we include also electron-electron scattering in the 
form of a mean-field potential which we determine by solving Poisson's equation 

d 2 V M F(z) e 

= -— [Pdopc(z) +Pe{z)\, (9) 
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with periodic boundary conditions VuF^z + d) = Vmf(z). e s is the absolute static background 
permittivity. Pdope(z) is the dopant density profile in the structure, and the electron density 
profile 

p e (z) = 2e£E -iG<^{E)r p {z)M*l (10) 

ap k 

with the factor of 2 from summation over the spin index. 

To determine the Green's functions, G ret (E) and G < (E), we loop iteratively between the 
quantum transport equations, Eqs. (f|) and (|5|), and the equations for the self-energies and 
mean-field potential Eqs. (||) - (|10|), until we reach a self-consistent solution for these equa- 
tions. To solve these equations, we assume that the system is an infinite periodic structure 
with period d. The boundary conditions G a p(E) = G a ipi{E + eSd) are applied between each 
period to model the bias drop per period in the structure, where a' and f3' are shifted one 
period to the left of a and (3. As starting conditions for the iterative loop, we assume the 
electrons are in a Fermi distribution with temperature T, and that the density of states in 
each Wannier subband is a simple stepfunction. However, once a converged solution was 
obtained for a given applied voltage, we used this solution as a starting point for the cal- 
culation at the next bias to facilitate the convergence process. To test for convergence, we 
compared the differences in the k-integrated Green's functions [G„* ,< (£ , )]j + i — [G^ c *' < (£')]j, 
and carrier density [n e ]i+i — [n e ]i, evaluated in two successive iteration steps, with a given 
tolerance value. A further test for convergence could be carried out by starting the calcula- 
tion at different bias points, e.g., starting at zero bias and increasing the voltage to a high 
bias point, then restarting the calculation at a high bias point and decreasing the voltage to 
a low value. 



C. Transport properties 

Having determined G ret (E) and G K (E), information about the stationary state of the 
system, such as excitation energies, energy renormalizations, broadenings or lifetimes, den- 
sity of states, populations, and distribution functions can be extracted. Examples of such 
information and their application are given in Sec. [HI A|. We can also evaluate the current 



density using the rate of change of the position operator z 

J(t) is the average current density flowing in the z direction, in system volume V. Recalling 
the separation of the Hamiltonian H into two parts, H Q and H' SCSutt , we can determine two 
separate contributions to the current J(t). The first contribution comes from H Q , 

Jo(t) = ~([H ,z]) = ^Tr{[H 0l z}G<(t,t)}, (12) 



which we evaluate in the energy representation 

2e y^y^ 1 

a/3 k 



J o = fjEE J^[Ho,zUG< aM (E). (13) 
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The second contribution is from H' 



^s'cattCO - Jsca,tt(t) + JMF(t) 

= ~^([H*xkt + HMF,z]), (14) 

where we evaluate Jmf similarly to J Q , i.e., Ref. [T8| 

2c { dE 

JlA¥ = TT^EE J ~^~l H MF, z]af3Gp aM (E). (15) 
a/3 k 

The current contribution driven by H scatt is written in the energy representation as (see 
Appendix 0) 

^att = ^EE/ ^^(^^(fO + GR^fOs^Scs)]^ 

a/37 k 

- ^[^(E)G^ k (E) + ^(E)G< atk (E)). (16) 

The superscript notation (a) in, for example, the self-energy T,^^(E) indicates that we 
take only the part of the self-energy which depends on G aa (E). The factor of 2 in Eqs. (|T3|). 
(|T5|), and ( ]T6] ) is from the spin index summation. 

We note here, that although we have divided the current into separate contributions 
driven by H a , Hmf, and i/scatt, this should not be taken to mean that the scattering processes 
play no role in driving, for example, the current J Q . Although, H sca tt does not appear 
explicitly in Eq. (0), it is instrumental in determining G < , and, hence, implicitly influences 
the current J Q . The scattering processes provide channels for energy dissipation from the 
electronic system, so that as the electrons descend downwards through the electric potential 
their kinetic energy does not increase but remains constant, and a steady current flow is 
maintained through the structure. 



D. Gain and absorption spectra 

We can also use the information contained in G ret (£') and G < (E) as a starting point to 
evaluate the gain or absorption spectra of the structure. To do this, we consider the linear 
response of the stationary state described by G vet (E) and G < (£') to a small perturbation 
from an optical fielcL From this response we determine the susceptibility x{ u )i an d from 
Maxwell's equations^ we obtain the gain coefficient 

^ _.Im[xH] i 
c n B 

where % ~ ^[T s is the background refractive index. The applied optical field 

E(r,t) = e z [ — S(oo)e lk ^ )y ~ luJ \ (18) 
J 2n 

which propagates in the y direction is included as a small time-dependent perturbation to 
the Hamiltonian: H — > H + SV(t), where the perturbation 

SV(r,t) = __l-^.A(r,t) + A(r,t)-p] + e0(r,t), (19) 
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with the momentum operator p, and m e (z) is the spatially dependent effective mass. The 
vector potential A(r, t) and scalar potential 0(r, t) are related to the optical field through 



J 2tc iuj 

4>(r,t) = 0, (20) 

in the Coulomb gauge. A discussion of the use of different gauges in this problem is given 
in Ref. For the results presented here, we apply the long-wavelength approximation, in 
which we neglect the spatial variation of the optical field across the structure. Hence, in 
this gauge, the perturbation is written in the energy representation as 



.e£(uj) 



to 



Pz 



m e (z) 



eS(uj) 



8V aP {u) = = —^-l[H ,zU. (21) 



a/3 



The linear changes in the Green's functions and self-energies caused by the additional 
term SV(t) in the Hamiltonian represent the linear response of the nonequilibrium stationary 
state to the applied optical field. These changes may be written asclJ 

5G vct (cu,E) = G Iet {E + tiw)[SV(u) + <5£ rct (cu, E)]G IC \E) (22) 

5G* dv {uj,E) = G adv (£ + hw)[5V(u) + 5£ adv (^ £)]G adv (£), (23) 

and 

5G < {u,E) = G TCt {E + hu) 5V{u) G < (E) + G<(E + %u) 5V(uj) G adv (E) 

+ G TCt (E + hu) 5E rc V, E) G<{E) + G rct {E + hu) 5£<(w, E) G adv (£) 
+ G<(£ + M<5£ adv (^£)G adv (£). (24) 

Note that 5G adv (u,E) ^ [5G iet (uj, E)}1 The expressions for 8E(E) take the same form 
as in Eqs. (|]) - (||) where the functions G(E) are replaced by 8G(E). For SH adv (E), the 
expression for <5£ ret (i?) is used with the replacement 5G < (E) — > —5G < (E) and SG ret (E) — > 
5G adv (E). 5'E < ^ TCt (E) must be evaluated self-consistently with 5G < / TCt (E), and therefore 
another iterative loop must be carried out, similar to the earlier procedure to determine 
G < / TCt (E) and H < ^ et {E). From the change in the Green's functions and self-energies, and 
using Eqs. (|TTD - (|16|), we obtain the change in current density 

5J(u) = SJ (u) + SJ MF (u) + SJ scatt (uj), (25) 

wherell 

Og /> cLE 

6Jo{u) = hV^^J ^ 0l zU SG^ k (u,E) + [5V,zU G^ k (u,E), (26) 



2g r- ^J^J 

6Jmf(u) = EE J -^[HmfJU 6G^ k (u>,E), (27) 

a/3 k " 
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and 



2e r dE ^ 

a/37 k 



G%^E + M 5Z%F>(u, E) + G%^E + hu) SE™(u, E) 



z i0 



-.ret (/J), 
J 77,k 



a/37 k 



, (28) 



where the factor 2 is again from the spin index summation. From the change in current 
density 5J(lj), we obtain the polarization 5P{uj) = iSJ(u)/u induced in the material, and 
hence the complex susceptibility 



6P(u) 



i 5J{ui) 



and from Eq. (|T?D we obtain the gain coefficient g{oS). 



(29) 



III. APPLICATION TO QUANTUM CASCADE STRUCTURES 

We have applied the theory presented above to some example QCL structures reported 
in the literature. These structures were grown in the GaAs/ALjGai-^As material system, 
but this theory is also applicable to other material systems. The examples we consider are 
[A] the first GaAs/Al^Gai-^As QCL structured reported by Sirtori et ai, with x = 0.33, 
and [B] a structured reported by Page et ai. with a similar layer design to structure A but 
with a higher AI content in the barrier, x = 0.45, resulting in higher barriers. This QCL 
structure operates in pulsed mode at room temperature. The lasing transition in structures 
A and B occurs between quantum well subbands. Other parameters for these structures are 
summarized in Table [TT]. In this section, we present and discuss results obtained for these 
structures. 



A. Nonequilibrium stationary state 

In this subsection, we describe some of the information about the nonequilibrium station- 
ary state of a QC structure under an applied bias that can be extracted from the Green's 
functions. In Fig. ^|, we show some examples of the diagonal elements of Im[G rot (i?)] for 
k = plotted as a function of the energy parameter E. Figure [|(a) shows two examples of 
this function evaluated in the Wannier basis. These examples were evaluated for structure 
A at a bias of 0.2 V/period. Each curve represents a Wannier state in one period of the QC 
structure. The Wannier states are not energy eigenstates of the system, and if a Wannier 
state is expressed as a linear superposition of the energy eigenstates, the position of the 
peaks in, for instance, the curve labeled u — 3 represent the energies of the eigenstates 
which comprise this superposition. 
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It is possible also to express the Green's functions in the Wannier-Stark basis. From the 



isis we obtain the trans- 
I Applying this transfor- 



diagonalization of the Hamiltonian H Q expressed in the Wannier 
formation matrix between the Wannier and Wannier-Stark basis.t 
mation matrix to G Tet (E) we obtain the example curves shown in Fig. 0(b). Each of these 
curves (diagonal elements of Im[G ret (i?)] for k = in the Wannier-Stark basis) represents 
a Wannier-Stark state. In contrast to the curves shown in Fig. 0(a) with their complicated 
structure of peaks, each curve in the Wannier-Stark basis consists of a single, broadened 
peak. This indicates that at this bias, the Wannier-Stark states approximate closely the true 
energy eigenstates of the system. The positions of these peaks Ei are shifted downwards by 
around 10 to 15 meV in comparison to the eigenenergies obtained from the diagonalization 
of H Q . This shift represents the renormalization of the Wannier-Stark energy levels due to 
the scattering processes described in H' Bcatt . The broadening of the peaks also originates 
from these scattering processes, and the width of the peaks gives a measure of the decay 
rate Tj of these levels due to scattering. In particular, we take Tj = A_E fwhm , where AE {whra 
is the full width at half maximum of the single peak in the function Im[G^ t k=0 (i?)], in the 
Wannier-Stark basis. In Sec. [111 Cj , we will use this information to analyze the gain spectra 
obtained from the NGF theory. 

The curves in Fig. were evaluated for k = 0. If instead we integrate over the in- 
plane wave vector k, we obtain the curves shown in Fig. 0(a). These staircaselike curves 
are the density of states (DOS) in two example Wannier subbands (the two Wannier states 
considered in Fig. 0(a) are the k = states in these subbands). The onset of each step in the 
staircase corresponds to the position of each peak in the corresponding curve in Fig. 0(a). 
If we carried out the k integral in the Wannier-Stark basis, we would obtain a single step 
instead of a staircaselike structure, corresponding to the single peak (for a given curve) seen 
in Fig. 0(b). If we sum the DOS curves for all Wannier subbands within a period, we obtain 
Fig. 0(b), which is the total DOS per period of the structure. The result in Fig. 0(b) is 
basis independent if we sum over contributions from all basis states within the energy range 
shown. 

We consider now the lesser correlation function G < (E). At equal times, G^ uk (t) = 

i(a£ k (t)a„k(£)) = i(h u \f:(t)), where (n^t)) is the occupation of a state k in subband v. If 
we move into the energy representation, and sum over the contributions from all states in 
all subbands in a period, we obtain the population n(E) = 2J2uk ^ uk (E) as a function of 
energy within the period. Figure 0(a) shows the normalised population per period per unit 
energy for several applied voltages. Dividing these curves by the total DOS curves for the 
corresponding bias gives the distribution function as a function of energy, i.e., 

t[h) ~ 2E, k Im[G- k (£)]- (30) 

The distribution function at 0.26 V/period in structure A is shown in Fig. f|(b). In this 
curve, the jagged points between and 0.2 eV are LO-phonon replicas, and the onset of a 
small population inversion at around 0.2 eV is also seen. 



B. Transport 

We consider now the transport properties of the QC structures described above. Figure |5| 
shows a plot of current density against V/period, calculated from Eqs. ( pi]) - fllCf). The 
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current density for structure B is reduced below structure A because of the higher barrier 



heights in B. The separate current contributions J Q and J sca tt + ^mf from Eqs. (]13"D - ( IB"]) 
are also shown. In both the structures, the main contribution is from J sca tt (<7mf is negligible 
in comparison with J sca tt)- This finding is in contrast to the behavior in simple superlattice 
structures where the dominant contribution to the current is from J . This difference could 
perhaps be explained by the far greater number of subbands within one period of a QC 
structure (in comparison to one period of a simple symmetric superlattice structure), in 
which the envelope functions ip a ( z ) are spatially displaced across the period. This opens up 
the possibility for many more scattering transitions which facilitates the transport of carriers 
across the structure through scattering processes. Comparing A and B, we see that J Q is 
almost identical for both structures. The current J sca tt in B is around half that in A so that 
the reduction in total current in B compared to A arises from a decrease in the transport 
driven by the scattering Hamiltonian H scau tt- 

To compare the theoretical results with experimental data reported in the literature, we 
show in Fig. || the results for A and B plotted against similar scales as the experimental 
plots shown in Fig. 3 in Ref. 21], and Fig. 4 in Ref. Comparing the results for structure 
A with the data given in Ref. 21, we see the voltage trend in both curves agree well. In 
the theory, the voltage tends to be lower for a given current than in the experiment, but 
this may be accounted for by an additional series resistance in the cladding layers. For 
structure B, we see again that the voltage trends agree well with the data shown in Ref. |22. 
In particular, when comparing the trends at different temperatures, we find that in both the 
theory and experiment the current is higher for a given voltage at the higher temperature. In 
both theory and experiment, the curves at both temperatures converge at a certain current 
density. Unlike the experiment, the two curves do not cross at the point of convergence in 
the theory, but diverge beyond this point. The theoretical curves reproduce the sharp break 
in the voltage trend at around 23 kA/cm 2 as seen in the experiment, and the theory shows 
that this break corresponds to the onset of a region of negative differential resistivity (NDR). 
NDR is not seen explicitly in the experimental data in Ref. [22] because of an artefact of the 
measurement systemic 

The break in the voltage trend is attributed^! by Page et al. to the breaking of the level 
alignment in the injector region and the active region causing an interruption in the current 
flow. To verify this, we show the wave functions and alignment of the injector level i and 
upper laser level u (in the Wannier-Stark basisil) at 77 K, for different applied voltages, in 
Fig. [7] Figures |7|(a) - (c) correspond to the voltages marked a, b, and c in Fig. ||(b). Figure 
|6](b) shows that as the voltage increases from 10 V, the current- volt age characteristic passes 
through a peak at 13 V, and then beyond this point the current decreases to a minimum 
at around 16.5 V. Comparing this behavior with Fig. |7|, we see that the change in voltage, 
and the behavior of the I-V characteristic, is accompanied by a shift in the position of level 
i from below to above level u. Level i passes through the resonance with level u which 
causes the peak in the I-V characteristic seen at 13 V. The population percentage in these 
levels is also shown in Fig. [7], and we see that away from the resonance position, ~ 50% of 
the population is in the injector level, with only ~ 10% in the upper laser level. Close to 
resonance, however, we see that the population is more evenly distributed between the two 
levels, with ~ 30% in the injector level, and 20% in the upper laser level. The wave functions 
in the near-resonance case (b) are less well-localized than in the off-resonance cases (a) and 
(c), and tend to spread across both the injector and active region. This increases the overlap 
of the two wave functions and facilitates the transfer of carriers between the injector and 
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active region. 

In the experimental data at 233 K, the break in the voltage trend is less pronounced, and 
this was attributed^ to the presence of electrons with higher energy leaking into continuum 
states, and maintaining the current flow. In the theoretical results, shown in Fig. |6](b), there 
is a small increase in current at 233 K, but it is difficult to make a quantitative comparison 
with the experimental data because of the uncertainty in the experimental data in the NDR 
regional Figure ^ shows the electron distribution functions f(E) at 77 K and 233 K, in 
one period of the structure at 14.4 V, i.e., just beyond the peak in the I-V characteristic in 
Fig. ^|(b). The thin horizontal line marks the conduction-band offset or barrier height. The 
peak in the distributions at around 0.2 eV corresponds to the population inversion in the 
upper laser level. The population inversion at 233 K is about a third less than that at 77 
K. In the high-energy tails, the distribution is slightly larger at 233 K than at 77 K. The 
difference in the distribution functions is not large, however, and although this difference 
may result in more electron leakage into the continuum at 233 K, it is not clear at this point 
if this is a sufficiently large effect to elucidate the experimental results. 



C. Gain and absorption spectra 

In this section, we apply the theory described in Sec. [II D| to structure A. We distinguish 
here between the gain coefficient g(u) evaluated using 5J{uj) from Eq. (p5|) , and the gain 
coefficient g {uS) which is evaluated neglecting 5J sca _ tt {uj) in Eq. (|25|) , and neglecting terms 
containing (5£</ ret / adv in 5G < (lo ) E) [Eq. (|24|)]. Hence, g {u) does not depend on changes in 
the self-energies, and it is simpler to evaluate since its evaluation does not require a further 
self-consistent calculation. 

In Fig. |9], we show the gain coefficient g{oj) calculated with the NGF theory for different 
applied voltages, between to 0.2 V/period. With zero applied bias, there is a broad 
absorption ranging from around 120 to 140 meV. As the applied voltage increases, the 
absorption gradually decreases in this range. There is a corresponding slow increase in 
absorption between around 80 to 100 meV. At around 0.18 V/period, a positive gain begins 
to appear in the range 120 to 140 meV, and this gain increases further as the voltage is further 
increased. These results were obtained in the Wannier basis and as stated in Sec. Ill A|, there 



is a problem in interpreting the origin of the different features in the spectra because these 
cannot be related to transitions between Wannier states. Later in this section, we discuss a 
way round this problem by considering a transformation to Wannier-Stark states. 

Before this we compare the results obtained for g{uj) from the more rigorous NGF theory, 
with the gain g (u) which is calculated neglecting terms from 5J scat t. This comparison is 
shown in Fig. [TT]. We see that the simpler theory gives rise to gain curves which are offset to 
higher values (i.e., absorption is reduced and gain is increased) in comparison to the more 
rigorous theory. 

To determine the origin of the different features in the spectra, we consider a different 
approach for evaluating the gain coefficient. In this approach, we transform the Green's func- 
tions G < / TCt obtained in the Wannier basis to the Wannier-Stark basis (G^/g Ct ) as described 
earlier in Sec. [Ill A|. From the diagonal elements of G^ s u , we obtain the populations of 



the Wannier-Stark levels, and from the positions and widths of the single peaks in Im[G^g u ], 
we obtain the Wannier-Stark level energies, Ei, for k = 0, and the corresponding decay rates, 
Tj. We treat each pair of Wannier-Stark subbands in a simple two-band model, in which 
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the optical response due to subbands i and j is given by,H 

lm[xij(uj)] = ^~ Y1 \ d ij\ 2 (Ak ~ /i,k) <$(^ + £ i)k - % k ), (31) 



e V 



k 



where d%j = ej dzip*(z) zif)j(z) is the dipole matrix element between the Wannier-Stark 
functions ipi(z) and if)j(z), and /^k and /^k are the nonequilibrium distribution functions 
in subbands i and j. The factor 2 comes from summing over the spin index. Assuming 
parabolic bands and replacing the 5 function with a Lorentzian to model the broadening of 
the levels we obtain 

Im[x(o;)] = — d Yl \dij\ 2 (ni-n>j)£ij(u), (32) 

(s/>Bi) 

with the Lorentzian C%j{u)) = {Yij/2ii) /[{fiu — AEji) 2 + (Tij/2) 2 ], where = Tj + Tj, 
and AEji = Ej — Ei (these parameters are defined in Sec. |111 A| ). rii (rij) is the 2D carrier 
density in subband i (j). In Eq. (|3"2"D, the contributions from all transitions with Ej > E^ are 
included. Using Eq. fl3"2p with Eq. (|I7|) gives us an estimate of the gain coefficient g ws {uj) 
within the Wannier-Stark picture. 



Disadvantages of this simpler approach, compared to the approach described in Sec. |IID 
are, firstly, that the Wannier-Stark levels become increasingly delocalized at low bias and 
an increasingly large number of basis states must be used to obtain a good approximation 
to the energy eigenstates as the applied voltage decreases. Thus, this approach becomes 
impracticable to apply at low bias. In addition, the function Im[G rct (i?)] is not always a 
simple Lorentzian as seen in Fig. ^|(b) for each Wannier-Stark level, but may, for some levels, 
have additional peaks or shoulders indicating that these levels do not approximate well a 
single eigenstate of the system but are superpositions of the eigenstates. This behavior 
occurs mainly at low bias, but can also occur when two Wannier-Stark levels are very close 
in energy which gives a double-peaked structure to lm[G ret (E)}, e.g., for the two levels in 
panel (b) in Fig. 0. In these cases, it is not possible to use this simple model to estimate the 
gain. Another drawback of this model is that by only making use of the diagonal elements of 
the Green's functions, G ws u and G^s w we have discarded information on quantum effects 
that are contained in the offdiagonal elements. 

In Fig. [11], we compare the gain coefficient g(u) from the NGF theory with the gain 
coefficient gwsi^) from the simple model in the Wannier-Stark basis. Spectra are shown 
for a bias of 0.12 - 0.2 V/period. We find that the simpler model reproduces well the main 
features of the NGF theory. As in the previous figure, we see that for the simpler theory 
in Fig. [II], the gain curves are again offset to higher values compared to the more rigorous 
theory. 

We consider now the origin of the different features in the spectra. This can be most 
simply understood if we look at the gain spectra obtained in the Wannier-Stark basis. For 



clarity, only these spectra are shown again in Fig. 12], and for a wider bias range, 0.1 to 



0.22 V/period, than seen in Fig. |TI| These spectra contain transitions between all possible 
pairs of levels with energy separation AEji i n the energy range shown on the x axis. At the 
lowest bias shown, i.e., 0.1 V/period, we see that the structure is almost transparent with 
only a very small absorption for the energy range shown. As the applied voltage increases, 
absorption increases in the range 0.08 to 0.11 eV. Simultaneously, gain appears and increases 
in the range 0.11 to 0.16 eV. The gain peak shifts to higher energies as the voltage increases. 
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To explain these features, we consider the curve at 0.2 V/period in Fig. |T2[ In Fig. [13|, 
we show the main Wannier-Stark energy levels involved in the transitions which give rise 
to the absorption and gain features seen at this bias. From an inspection of the transitions 
contributing to the sum in Eq. (p2[), we find as expected that the gain feature originates 
mainly from the transition (marked B in the figure) between the upper laser level (labeled 
2) and the lower laser level (1). There is also some contribution to the gain from transitions 
between an injector level (2') and the lower laser level (1), and between other injector levels 
and the upper or lower laser level (not shown). The absorption feature below 0.1 eV is due 
mainly to transitions (marked A in the figure) from the upper laser level (2) to levels in the 
continuum (3 and 3'), with additional contributions also from transitions between injector 
levels. The relative population in each level is also shown in the figure, and levels 2 and 
2' have a greater fraction of the population compared to levels 1, 3, and 3'. This supports 
the attribution of gain to transition B, and absorption to transition A. The largest fraction 
of the population ~ 70% is in the lowest injector level (not shown in the figure). Level 
contains around 2% of the population. 

When we consider the relative populations in these levels for different applied voltages 
(in the range 0.1 - 0.22 V/period as shown in Fig. |l|), we find that the population in the 
lower laser level remains at around 1% at all voltages. The population in the upper laser 
level, however, changes from around 0.2% at 0.1 V/period to over 10% at 0.22 V/period. 
Hence, there is a transition to a population inversion between the two levels as the bias 
increases. The population in level ranges from around 20% at 0.1 V/period to around 2 to 
3% above 0.16 V/period. The largest proportion of the population remains in the injector 
region, particularly in the lowest injector level which forms a reservoir of around 60 to 70% 
of the carriers at all applied voltages. 

We also note here that the broad gain or absorption features seen in Fig. [12| are not due to 
transitions between a single pair of levels but contain contributions from many transitions. 
For example, in Fig. |TJ], we show the different Lorentzian contributions (thin lines) to the 
gain curve (thick line) at 0.22 V/period. Although the main contribution to the absorption 
feature (thin line labeled A) originates from the upper laser level - continuum transition, 
and the main contribution to the gain feature (thin line labeled B) arises from the upper - 
lower laser level transition, there is also a substantial contribution from other (energetically) 
neighboring transitions. The figure also shows that the inhomogeneous broadening due to 
the contributions from different transitions, and the broadening due to scattering, observed 
in the linewidth of each individual Lorentzian, are comparable in size. At higher applied 
bias the number of contributing transitions decreases, e.g, at 0.3 V/period, the contributions 
to the gain feature come mainly from only two transitions. 

Finally, in the last set of results, we consider structure B. Experimental data showing the 
light output vs. current density is shown in Fig. 4 of Ref. At 77 K, the light output is 
seen to increase with increasing current density until it reaches a maximum at around 22 
kAcm -2 , corresponding to the peak in the I-V characteristic seen in Fig. §|(b). The light 
output goes to zero for higher current densities. A similar behavior is seen at 233 K but with 
a much reduced light output. Figure [TJ shows the gain coefficient g ws {^) calculated at 77 K 
for structure B at different applied voltages around the peak of the I-V characteristic shown 
in Fig. §|(b). The gain feature increases as the applied voltage or current density increases, 
reaching a maximum around the peak of the I-V characteristic, and then decreases beyond 
this point. A similar behavior is also seen at 233 K except that the gain is reduced at 
the higher temperature. For comparison, a gain curve calculated at 233 K, at the peak 
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current density 22 kAcm -2 , is also shown in Fig. |15| The gain feature in this curve is much 
smaller than the gain feature in the corresponding curve (at 22 kAcm -2 ) at 77 K. Hence, 
the behavior of the calculated gain curves correlate well with the light output vs. current 
density curves measured in the experiment. 



IV. CONCLUSION 

We have applied an NGF theory to obtain a description of the nonequilibrium stationary 
state of QC structures under an applied bias. Using this information, we evaluate the 
current- volt age characteristic of example QC structures reported in the literature. The 
theoretical results are quantitatively close to experimental I-V data, and reproduce well the 
trends seen in the data. In addition, we determine two contributions to the current density. 
The first contribution J Q is driven by H Q which is the Hamiltonian for the superlattice 
potential with applied bias. The other contribution J s ' catt is driven by H sca ,tt + ^mf which 
describes the scattering processes in the structure. We find that, in contrast to simple 
superlattice structures, J s ' catt is the main contribution to the current in the QC structures 
we consider. 

In addition, we have extended the theory to determine the linear response of the nonequi- 
librium stationary state of these structures to a small applied optical perturbation. This 
enables us to evaluate the linear susceptibility and hence the gain or absorption spectra of 
the structure. We compare the spectra obtained using a more rigorous NGF theory in which 
the changes 5G, 5X1, and 5J scatt due to the optical perturbation are considered, to simpler 
models in which (i) only 5G is considered, or (ii) by summing over transitions in a simple 
two-band model (summing over different pairs of bands) with Lorentzian broadened levels. 
We find that the simpler models result in spectra which are offset to higher values than in 
the more rigorous theory. We have also made a detailed analysis of the origin of the different 
gain and absorption features in the spectra, and of the redistribution of population within 
the Wannier-Stark levels as the applied voltage changes. The gain and absorption features 
correlate well with the distribution of population within these levels, which is determined 
from the NGF theory. 
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APPENDIX A 



The construction of the Wannier function basis used in the calculations is described in 
this appendix. As the first step, we solve the one-dimensional Schrodinger equation 



h 2 d 2 

2m e (z) dz 2 



+ V(z) 



4>(z) = Eif>(z), 



(Al) 



for the envelope functions ip( z )- The spatially dependent superlattice potential V(z) and 
the effective mass m e (z) are assumed constant in each semiconductor layer, i.e., V(z) = V 
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and m e (z) = in the barriers, and V(z) = and m e (z) = m w in the wells. Equation (|A1|) 



can be solved with a transfer matrix method (for a textbook discussion see, e.g., Ref. pop . 
In this approach, the envelope function in a semiconductor layer j is written as 

^j(z) = Aje^^-^ + Bje- lk ^ E){z -^\ (A2) 



where Zj labels the position of interface j, and kj(E) = y2rrij(E — Vj)/h with rrij and Vj 
the mass and potential in that layer. Applying continuity conditions 

tpj = tpj+i, (A3) 

1 dVj _ 1 / A4 n 



at the interface j + 1, gives 



with 



rrij dz 




M j (E)[J , (A5) 



jV ' 2 Wl - 2z±l*i , ) e ifcj(^+i-*j) fi + I2i±i*n e -iM*i+i-*i) /' 1 ; 

If a single period of the structure contains M semiconductor layers, the Bloch condition 
+ d) = e iqd ip q (z) implies 

= TT MA FA ( 



■j+i fc j \ c —ikj(zj+-\ —Zj) 




Y[M J (E)(^)=e^[^\. (AT) 




3=1 

For a given value of q, only certain values of E allow the solution of Eq. (|A7j), and this 
defines the miniband structure E u (q). For each q value, we determine E u (q) numerically by 
looking for the zeroes of the determinant of the matrix product Y\jL 1 Aij(E). These zeroes 
were found by stepping through the energy E, and comparing the signs of the determinant 
for two consecutive values of E (separated by AE). When these signs are opposite this 
sets the first coarse bracketing of the zero position. This position is then further refined by 
halving this interval and successive intervals up to a hundred times, while comparing the 
determinant signs of the interval endpoints at each iteration. Once E u (q) is determined, we 
can also obtain A u (q) = A\ and B u (q) = Bi, and the Bloch functions ipg(z), from Eqs. ( |A2[ ), 



and (|A7[ ). In the calculations reported here, E u (q) was evaluated for 500 q points for 



each miniband, and AE was set to 5 meV. The Bloch functions were evaluated on a position 
grid with 1500 points per period. In structure A (B), there are eight (nine) minibands below 
the conduction-band offset. For the calculations here, we included one miniband above the 
continuum, so the band index v runs from 1 to 9 (10) for structure A (B). 

In the next step, we construct the Wannier functions (associated with miniband v) 

tfj^(z - nd) = \ — dq e~ inqd (A8) 

V 2ll J—ir/d 

from a superposition of the Bloch functions in miniband v. The Wannier functions con- 
structed by this superposition are not unique, and can be very different depending on the 
phase of the Bloch functions which can be chosen arbitrarily for each value of q. Ideally, we 
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would like the Wannier functions to be as spatially localized as possible to reduce the in- 
teraction matrix elements between Wannier functions in different periods. As a first step to 
construct Wannier functions with this property, we first fix the phase of the Bloch functions 
such that these functions are real at some arbitrary position x s . Different values of x s can 
be tested, and the resulting Wannier functions checked to see if they are fairly well-localized, 
e.g., within one period of the structure. We then express the position operator z in the basis 
made up of this initial set of Wannier functions. Finally, we diagonalize the resulting matrix 
representation of z, and the resulting eigenfunctions give us the required set of Wannier 
functions associated with a given miniband. This process is repeated for each miniband. 

The Wannier states associated with a given miniband are degenerate, and their energy 
expectation values lie at the center of the miniband. Their wave functions are spatially 
displaced from each other, with a separation given by the period of the structure, i.e, there 
is one Wannier state in each period of the structure. To set up the matrix representation of 
z for a given miniband, we use Wannier states in 11 neighboring periods. The large number 
of periods used in this construction is necessary to improve the numerical accuracy of the 
result. 

In the transport and gain calculations, however, we consider only couplings and ma- 
trix elements between Wannier functions in the same period, and nearest-neighbor periods. 
Hence, the matrices representing the Green's functions and self-energies are constructed with 
27 basis states (for structure A) from three periods. 



APPENDIX B 



This appendix describes in more detail the evaluation of the self-energies and scattering 
matrix elements in Eqs. (Q) - (§). The self-energies contain summations over k' of the form 
iV^^k, h')\ 2 {G.F.} where V a p is a generic matrix element representing interface rough- 

k' 

ness, impurity, or phonon scattering, and {G.F.} represents a Green's function Gpp^i(E'), 
with E' = E or E' = E ± E phon , and including, for the case of phonon scattering, a phonon 
distribution factor. Taking the summation to the continuous limit leads to 



J2\V a p(k,k')\ 2 {G.F.} 
k' 



A poo r2n 

— / Q dk'k'l dO\V a p(k,k', 



\G.F.}. 



(Bl) 



We carry out the angle integration assuming the Green's function term does not depend on 

d9f(k',9) we obtain 

10 



the angle, and defining the angle-integrated quantity [f(k')]$ 

A 



ElK,(k,k')| 2 {G.F.} 
k' 



(27T) 2 

A m e 

1 Ap 
2tt 2 



dkk 



V^(fc,fc')r AG.F.} 



dEy 



V a p(k,k')\ 2 AG.F.} 



V a p(k typj k>)\ 2 ]ldE k , {G.F.} 



(B2) 



The integral over momentum k' has been transformed to an integral over energy E^, and the 
density of states p Q = m e /(iih 2 ). In the last line of Eq. (|B2|) we assume the matrix element 
V a p(k, k') is slowly varying compared to the Green's function term and can be taken out of 
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the integral. It is evaluated at fixed momenta k tJP and k' typ , and the choice of these momenta 
is described below. 

From a comparison of the factors outside the integral, in the last line of Eq. (|B2|), with 
a scattering or transition probability rate derived from Fermi's golden rule we can defmeU 
a fictitious scattering rate (in energy units) 

Ap r 

la/3 = 



2 



\V a p(k typ ,k[ yp )\ 2 \ e . (B3) 
With this definition we can rewrite Eqs. (|Bl| ) and (|B^) as 

£|y Q/3 (k,k')| 2 {G.F.} = M / d E k , {G.F.}. (B4) 

k , 27T J 

Estimates of 7 for different scattering processes are given in Table |. 

As stated after Eq. and earlier in this appendix, we evaluate the scattering matrix 
elements V$ oa/TOUsh/,mp (k,k') using fixed momenta k typ and k' typ (with the corresponding 
energies E typ and E' t ) to accelerate the numerical computation. To fix E typ , we consider the 
energy dependence of the scattering matrix elements. For LO-phonon scattering, there is an 
energy threshold for phonon emission because of energy conservation and the fixed phonon 
energy Ej^q. The scattering matrix element is maximum at the emission threshold, and 
decreases monotonically with increasing energy above this threshold. To obtain an estimate 
of the scattering matrix element, which lies between the higher values near threshold, and 
the lower values far above threshold, we set E typ one LO-phonon energy above the LO- 
phonon threshold. We then fix E' t = E typ + AE a/3 — E LO , where AE a/3 = [H ] aa — [H ]pp 
(see Fig. [L6|). To test the sensitivity of the results to the value chosen for E typ , we have 



carried out runs with other values of -Etyp, e.g., at threshold, |-Elo above threshold, and 
2E\ J o above threshold. The calculated current density changes by at most ~ 9% (at some 
bias points with E typ at threshold) but in most cases the change is around 5% or much less 
(1 - 2%). This tends to support the assumption that, for the LO-phonon process, the results 
are not very sensitive to the specific value of -Etyp- For impurity and interface roughness 
scattering, we set E typ = 15 meV, and E[ = E typ + AE a p. Test runs were also carried out 
for E typ = 1,7, and 30 meV. The results are more sensitive to the value of Etyp- For E typ = 7 
and 30 meV, the calculated current density changed by at most 15%. For E typ = 1 meV, 
near the bottom of the subband, the difference was much larger, ranging from 10 to 50%. 
The value E typ = 15 meV was chosen as a value that lies near the centre of the distributions 
in each subband, to give an estimate of the average scattering matrix element. 

With the momenta fixed, the scattering matrix elements can be taken outside the integrals 
in Eqs. @ - (||) , and the self-energies depend only on integrals of the form [as shown in Eq. 

poo 

Gf a {E) = / dE k ,G aa v(E), (B5) 
Jo 

and 

roc 

hn[G< a (E)} = / dE k> lm[G< a ^{E)]. (B6) 
Jo 

Note that the diagonal elements of G < are pure imaginary. A problem arises in evaluating 
Gaa{E) when a momentum- independent scattering matrix element is used because this leads 
to a divergence in the integral as E k > — > 00. To deal with this problem, a cutoff energy for 
the upper limit of the integral is used in the numerical integration. The following subsections 
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give more detail concerning the evaluation of the matrix elements for the different scattering 



processes (see, also, Refs. [27]. |28| ). 



1. Interface roughness and impurity scattering 

To derive matrix elements for interface roughness scattering we consider an interface j 
located at z = Zj with thickness fluctuations £j(r) of the order of one monolayer. We assume 
correlation functions for the fluctuations given by 

(fc(r)> = (B7) 
(&(r)fr(*0> = ^/(|r-r'|) with /» = r? 2 exp (-£) (B8) 

7] denotes the root-mean-square of the roughness height and A is a typical island size. Cor- 
relations between neighboring interfaces are neglected. The angle brackets [as stated after 
Eq. (K)] denote an average over different distributions of thickness fluctuations. 
The Hamiltonian for interface roughness scattering is written as 

4ou gh = E [K^f,(p)^(k + P)a n ,,(k) + H.c.], (B9) 
k, P 

mfi,ni/ 

with the matrix element 

V^tM =Y,jJd 2 r e-ni(r)A^(^ - md)Mzj ~ nd), (BIO) 

where AE C is the band offset. 

Within the self-consistent Born approximation, the self-energy contribution from interface 
roughness scattering is written as 

KaS(E) = E (C^ h (k - k')^° Q u f(k' - k))G<- k T h (E). (Bll) 

/3/3',k' 

This equation is more general than Eq. (^) since it includes the offdiagonal contributions. 
Now we assume that the diagonal parts of G^p°y dominate, and we keep only the terms 
(3 = (3' in the summation. Then we obtain for the matrix element term 

(K7 h (-p)^r h (p)) = &(-p)ct,(p)> 

= ^f^<E / d 2 r 1 e i P^T 1 )r» 1 (zf-m 1 d)^z f -nd) 

xE/ d 2 r 2 e-^ r ^ 3 {v 2 )r u {z ] -nd)ij^{z 3 -m 2 d)] 

j 

{AEc)2 J d 2 n J d 2 r 2 e^-^f(\ ri - r 2 |) 



A 2 

x E^i( z i ~~ mid)\ip v (zj - nd)\ 2 ip^{zj - m 2 d) 
j 

A2 - J d 2 re^f(\r\) 
x E^i( z i ~ rnid)\ip v {zj - nd^ip^Zj - m 2 d), (B12) 



(AE, 
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where we have expanded the general indices ax,ac2,{3 in terms of the period and Wannier 
level indices [see after Eq. (Q)], and defined p = k' — k. Because of the orthogonality of the 
wave functions, and because the wave functions extend over many interfaces, the sum over 
j tends to vanish for m x ^ m 2 and /ii 7^ fi 2 - This suggests replacing this term in the last 



line with F^_ n 5 mi ^ m2 5^^ 2 where 



(B13) 



Applying Eq. ( p8|) we obtain 

/yrough/ xyrough/ \\ 



(AE C ) 2 V 2 



poAEy 



E P /E X )W 



HE; 



pi- 



(B14) 



with E x = h 2 /2m\ 2 and E p = h 2 p 2 /2m = E k + E k > — 2\/E k E k > cos6, where 6 is the angle 
between k and k'. The subscript 1 from the indices mi and /ii is dropped for simplicity. 
Substituting this latter result in Eq. ( Pll|) gives Eq. (§) which contains only the diagonal 
terms in the self-energy and the Green's function. We now follow the procedure outlined 
in Eqs. ( Bl|) - (|B2"|) and take the summation over k' in Eq. (||) to the continuous limit. 
We observe that both the self-energy and the Green's function depend only on E k and E k >, 
but not on the angle 9. Hence, as shown in Eq. (|B~2"|), the angle integration over the matrix 
element can be carried out analytically, and we define the angle-integrated quantity (Ref. |29| , 
2.575): 



lK^f,(k-k')| 2 ) 



2k 



d6 T[E k + E, 



with 



1 + 



(AE C )V 
PoAEx 



E k + E k i 




E\ E x 
E(x) is the complete elliptic integral of the second kind which is of order 71/ 2 
E(x) > E(l) = 1 (Ref. |I| 17.3). Therefore we set E(x) w tt/2, and fixing E k = 



(B15) 



(B16) 

= E(0) > 
E typ and 



E k i = E typ + AE al 3 as described earlier, we define the fictitious interface roughness scattering 
rate [using Eq. ( p3|) 1 



rough 
I(m—n),ij,u 



ii (AE, 



1 9 



E\Fi^_, 



E x + 



E t yp - 



E typ + AE 



a/3 



E X + (/#typ 



E tyP + AE, 



(B17) 



a/3 



Impurity scattering is mediated through the (3D) Coulomb interaction V ~ 1/q 2 ~ 
l/|k— k'| 2 . Neglecting the angle-dependent term gives \V\ 2 ~ l/(E k +E' k ) 2 . Fixing^ = E typ 
and E k i = E typ + AE a p as before, and neglecting factors of 2, we obtain an order-of- 
magnitude estimate of the impurity scattering rate 

'jdz\P a {z)M*)\ 



imp 

la/3 



7- 



imp 



par (l + AE a0 /E 



(B18) 



typ) 



The parameter 7p™ p (see Table |) is estimated from calculations of impurity scattering in 
simple superlattice structures.Hl 
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2. LO-phonon 



To derive the electron-LO-phonon scattering matrix element we start from the interaction 
Hamiltonianii 

Hio = E ^u(Q)ie- ^Q ■ R tf Q - e^%], (B19) 
Q 



with 



«(Q) 



e 2 huj] n ( 1 1 



(B20) 



.2V Q 2 Veoo e s 

Q = (Q) 9z) is the 3D phonon wave vector with the 2D in-plane component q and the 
component q z in the growth direction. Similarly, the 3D position vector R = (r, z). 6q (6q) 
is the phonon creation (annihilation) operator, hu\ = E\ Q is the LO-phonon energy, and V 
is the crystal volume, and e s are the high-frequency and static absolute permittivities. 
The matrix element 

V% = H#i |/3) = Y,MQ)(^\e- iQ -%-e^%\(3), (B21) 

Q 

is written with initial and final states, |/3) and given by 

= |# ki/3 (r,z))|n(u; lo )), 
I a) = |*k',a(r,z)>|nW±l). 

The ket containing n(uJi Q ) is a phonon number state. The upper (lower) sign in the ket 
\n{uj\ ) ± 1) corresponds to phonon emission (absorption). The ket containing \l/ gives the 
electron state, and the electron wave function ^^(r, z) = A~^e % ^' r ^)fi{z). Evaluating Eq. 
( [B2TD leads to 

|V^(k,k')| 2 = \{a\KW = E« 2 (Q)[^K) + |±|]l^fe)| 2 5kMc Tq , (B22) 

Q 

with 

M aP {q z )= / dz e™* z i>* a (z)Mz)- 
Jo 

L w is the distance in the ^-direction over which the wave functions ip(k z , z) extend. Con- 
verting the sum over Q in Eq. (|B22|) to an integral, and evaluating the in-plane component 



q with the help of the Kronecker delta gives 

V^Xtf = [n M + i ± ^jy q J^01, (B23) 

where p = |k — k'| 2 = k 2 + k' 2 — 2kk' cos 9. 

As described earlier in this appendix, the matrix element is evaluated within a self-energy 
integral, for instance of the form: S(k, E) = J2w k')| 2 {G.F}. Following again the 

procedure shown in Eqs. (|BT| ) - ( p2|) , we take the summation over k' in the self-energy to 
the continuous limit, and we define the angle-integrated quantity 
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lOk,k')| 2 



/■2tt 

C / d# 



2nC 



C r°° 



dq 



Po 



dq z 



*> |M Q/3 fe)| 2 
oo 9z q 2 + k 2 + k' 2 -2kk' cos 9 
\M af3 (q z )\ 2 

(q 2 z + k 2 + k' 2 ) 2 - 4k 2 k' 2 
\M a p(q z )\ 2 



2m 



+ E k + E k ,) 2 - AE k E k , 



where C = E lo e 2 [n(u lo ) + \ ± ±]/(4vre p .4). 

Using this result in Eq. (p3|), and fixing the energies E k = 
described earlier, we can define the scattering rate 

io r / i , ii^ioe 2 /°° , \M af3 (q z )\ 2 



E typ and 



4-Etyp-Etyp 



(B24) 



E 'ty P as 



(B25) 



3. Acoustic phonon 



The phonons are implemented as an artificial optical phonon with a phonon energy E ac 
which should be smaller than /c^T and which should not be commensurable with the optical 
phonon energy. The matrix element, or equivalently, the fictitious scattering rate 7^, is set 
to 



la/3 



par 



dz\ip* a (z)ip/3(z)\ 



(B26) 



APPENDIX C 

A derivation of the current contribution J sca tt [Eq- QIBD 1 is given here. We assume that 
the scattering Hamiltonian H Bca tt has the form 

O a ^{t)d\^ s {t)d^ s {t), (CI) 

a/3 

where O a p{t) may be just a scalar time-independent matrix element, i.e. O ak „3k'(0 = Vak^k', 
as in interface roughness or impurity scattering. Alternatively, as in phonon scattering, it 
may be a time-dependent operator containing the phonon operators b{t) and W(t). Then we 
find 

^scatt(^) = 77 7" ( [-^scatt j ^] ) 

v n 

= ^7 E i («ik S WP«k, 7 k'(^ 7 ^ - z ai d lK ^{t)]d^ s {t)) . (C2) 

k,k',s 

In the following, this expression will be expressed in terms of Green's functions and self- 
energies. For this purpose, we consider the contour-ordered Green's function: 

F c aM (ri,r 2 ) = -iT c {(a / 3 k ' s (ri)a^ ks (r2)6 Q k, 7 k'( r 2)2;7/3 - ^ 7 7 k, / 3k'( r i)a/3k's(ri)a) iks (r2))}, 

k,k',s 

(C3) 
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where T c orders the times arguments r on the complex contour from r = — oo + i0 + to 
r = — oo — i0 + (see Ref. [IT]). By taking the lesser component Fl (37 (t,t), one obtains the 

k,k',s 

summand in Eq. flC2|) . 

Using the interaction picture with superscript D, F c is now evaluated according to the 
standard perturbation expansion of 



F%y (r^Ta) = -ifJ(expi [ --jdTH D (T)j a^ s (ri)a^ s (r 2 )0 Q k,7k'( r 2)-27/3 



- 2a70 7 k^k'(n)a^k' S ( r l) a ak S ( T 2 



(C4) 



The exponential term can be expanded as 



exp [-- I ,hH L '(T\ 



dr 0«5 P ,eq(r)af p t s ,(r)a^(r). 



(C5) 



Substituting this expansion into Eq. (|C4|), and noting that terms containing only a single O 
are zero after averaging, the lowest-order nonvanishing terms give: 

F c a01 (ti,t 2 ) « - J] / rfr [G^^n^O^^ 



k.k'. 



with the bare Green's functions G^ ks g ks 



7 k ) /3k'(Ti)G^ k / S)5k / s (r 1 , r)0 5kVk (r)G^ ks aks (r, r 2 ) (C6) 



ri,r) = -«T C { (03^(71)0^^))}. Only Green's 
functions diagonal in the momentum and spin indices are kept in Eq. ( |C6| ) . In order to be 
consistent with the perturbation expansion in the Green's functions, further terms are taken 
into account, which replace the bare Green's functions by the full Green's functions. Then 
we find 

k,k' k ' k '' s n 5,k J 

- Za^ySk ( T l' T )^5ks,aks( r 5 T 2)] (C7) 



<5,k 



where 



ek' 



denotes the part of the self-energy which exhibits O ai on the right-hand side, and 
E ?a?k ( r i> r ) = Z^ 7k,/3k'(ri)G^ k , Si£k , s (ri,r)O ek / i5 k(T-) 

ek' 



(C8) 



(C9) 



with 7 /3 on the left-hand side. In the second term in Eq. ( |C7| ) we have exchanged the 
dummy indices 6 and e, and in the first term we have exchanged k and k'. 

For diagonal self-energies, which depend on diagonal Green's functions (see Sec. 



one has £^ k 



Ss^—jit where £™ k is defined after Eq. flIS|) . Using Langreth rules 
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and changing to the energy representation leads to Eq. ([II]) which is given again here for 
reference: 



Jscatt = kE/^fe(<f( £ ) + Gt±{E)V<%{E))z lP 



&$(E)G%(E) + KT( E ) G ^(E)) . (CIO) 



We can obtain some insight into this expression for J sca tt if we consider only the diagonal 
terms of the above equation, i.e., we set 7 = (3 in part I, and 7 = a in part II, to obtain 

/ dE I — fdE( Gf p s; d / (Q) + G% ) zpp, (Gil) 

la lb 

JdEll J dEz aa (E<®G^ + S r ^G< Q ), (C12) 



iia 



where for brevity we neglect the index k, and £ and G are understood to be functions of 
E. We observe that each term in the integrals above is a product of a Green's function 
and a self-energy, i.e., of the form G£ or EG. These terms can be interpreted as scattering 
rates. In particular, Eq. flCH| ) contains information about scattering rates into and out of 



the state (3, and Eq. ( |U12| ) describes scattering into and out of the state a. To be more 
specific, we interpret the terms containing E ret or S adv as a scattering-out rate, e.g, the term 
la, / d£GlE^ v ' a ', can be interpreted as a rate Tg^ a for scattering out of state (3 into 
state a. Similarly, term lib, / dEY7*^G< a , is interpreted as a rate r°^a for scattering 
out of state a into state (3. On the other hand, the terms containing S < are interpreted as 
scattering-in rates. Thus, term lb describes the rate rj^Cg for scattering into state (3 from 
a, and term Iia describes the rate FgV^ for scattering into state a from (3. We note here 
that Y°^t a = r'f_+ a < and P a n _^ = r° u _^ > 0. Combining Eqs. (|CTID and (Pl2l) gives 

JdEi + n = (r°£ a + - z aa (r^ a + r^y 

= (r^ + CiT)(^-^). (cis) 

Thus, this expression is the product of the distance \zpp — z aa \ with the net transfer rate 

(Pp^la + r^i^f) between state (3 and a, and we can interpret this as a velocity or charge 
transfer rate from e.g., zpp to z aa , i.e., a current flow from zpp to z aa . (The direction of 
current flow depends on the net transfer rate.) 
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Interface roughness 


LO-phonon 


Acoustic phonon 




& impurity 






intrasubband 


0.05 ps 


0.2 ps 


> 7 ps 


intersubband 


> 0.2 ps 


> 0.6 ps 


» 8 ps 



TABLE I: This table gives order-of-magnitude estimates of inverse scattering rates between Wan- 
nier levels in structure A (see Sec. Ill) at 77 K. The times given in this table are obtained from 
h/j with 7 (in energy units) defined in Eq. flB3|) . More specifically, for interface roughness and 
impurity scattering, Y^-n) ^u^^ap f rom Eqs. ( |B17j ) and QB18 ) was used, and for acoustic phonon 
scattering, Eq. ( |B26| ) is used. For LO-phonon scattering, Eq. ( |B25| ) without the phonon distribu- 
tion factor was used. The parameters used in obtaining these estimates are as follows. For interface 
roughness scattering: AE C = 0.27 eV, r] = 0.28 nm, A = 10 nm. For impurity scattering: 7p™ r p = 5 
meV. For LO-phonon scattering: E\ Q = 36.7 meV, e s = 13.1eo, £oo = 10.9eo- For acoustic phonon 
scattering: 7^. = 0.1 meV. 
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Structure 


X 


n e (cm 2 ) 


d (nm) 


N p 


A (Ref. P) 


0.33 


3.9 x 10 11 


45.3 


30 


B (Ref. 


0.45 


3.8 x 10 11 


45 


36 



TABLE II: Parameters for example GaAs/AL E Gai_. r As QC structures: x is the aluminium content 
in the barrier, n e is the 2D carrier density in one period, d is the period length, and N p is the 
number of periods. 
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FIG. 1: Example of the conduction-band lineup in a quantum cascade structure with an applied 
bias. The arrows indicate the direction of electron flow in the structure. 
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Structure A 77 K 



0.2 V/period 




E(eV) 



FIG. 2: Examples of the diagonal elements Im[G^ k=0 (£')] in one period of structure A with an 
applied voltage of 0.2 V/period. a) Wannier basis, b) Wannier-Stark basis. 
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FIG. 3: a) Examples of density of states (per period) in Wannier subbands: ImfG^y (i?)] = 
2Ek lm [ G l%k( E )\- b ) Total density of states per period, £ I/ Im[G2*(.E7)]. 
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FIG. 4: a) n(E), normalized population per period for different applied voltages, b) f(E), 
distribution function in one period at 0.26 V/period. The thin vertical line marks the position of 
the conduction-band offset between the thin (first) well in the active region in this period, and the 
thick barrier separating this well from the preceding period. The bottom of this well is at around 
eV. The peak at around 0.2 eV corresponds to a population inversion in the upper laser level in 
the active region. 
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V/period 



FIG. 5: Current density vs. voltage/period for example QC structures. Structure A: Thick lines. 
Structure B: Thin lines with symbols. 
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FIG. 6: Current-voltage characteristic for comparison with experimental data for structures A and 
B (see text for references). The alignment of the injector level and upper laser level at the positions 
marked with arrows in (b) is shown in Fig. 0. 
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position (nm) 



FIG. 7: Wave functions (mod. squared) of the injector level i (solid line) and upper laser level 
u (dashed line) in structure B for different applied bias at 77 K. The energetic positions show the 
alignment of these levels. Parts (a) - (c) correspond to each of the voltages marked with arrows in 
Fig. |6](b). The percentage of population in each level is also shown. AE is the energy separation 
between the two levels. 
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0.002 0.004 0.006 

f(E), dist. func. 

FIG. 8: Electron distribution functions f(E) at 77 K and 233 K in one period of structure B at 
14.4 V. 
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FIG. 9: Gain and absorption spectra g{oj) for structure A with different applied voltages, calculated 
with the NGF theory. 
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FIG. 10: Comparison of g{oj) and g (w), with different applied bias, for structure A. a: solid lines, 
b: dashed lines, c: dotted lines, d: dot-dashed lines, e: dot-double-dashed lines. 
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FIG. 11: Comparison of g(to) and gws(w), with different applied bias, for structure A. a: dashed 
lines, b: dotted lines, c: solid lines, d: dot-dashed lines. 
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FIG. 12: Gain spectra ffws(w) evaluated in Wannier-Stark basis for structure A. 
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position (nm) 

FIG. 13: Wannier-Stark levels in structure A at 0.2 V/period. 
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^photon 

FIG. 14: Contributions of individual transitions (thin lines) to gain curve gws(w) (thick line) at 
0.22 V/period for structure A 
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FIG. 15: Gain curves gws(w) at 77 K for structure B, for different applied voltages. Thin solid 
line with symbols is at 233 K. 
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AE ap > E L0 > 




FIG. 16: Examples of -E'typ and E' typ selection, (a) LO-phonon scattering for the case AE a /3 > 
Elo > (requires E typ > 0). Similar figures can be drawn for the cases E^o > AE a p > 
(requires E'typ > -Elo — AE a p), and AE a p < (requires E typ > \AE a p\ + Elo)- For all three 
cases, E' typ = E typ + AE a p — Elo- (b) interface roughness and impurity scattering for AE a p > 0. 
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